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Preface 


Two  complementary  analyses  of  the  time-harmonic  scattering  by  a  penetrable  wedge 
are  presented.  The  distance  from  the  apex  (appropriately  scaled  by  the  wavenumber  in  the 
exterior  region)  of  the  exciting  line  source  is  the  single  length  scale  in  this  infinite-domain 
boundary  value  problem.  The  work  summarized  herein  represents  two  mathematical  ap¬ 
proaches  (among  a  series  of  candidates)  to  solve  this  important  scattering  problem  and  to 
visualize  the  wave  physics. 
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Figure  1:  Wedge  Geometry 

1  Introduction 

Consider  the  two-dimensional  scattering  of  a  time-harmonic  sound  wave 
generated  by  a  line  source  and  incident  upon  a  penetrable  wedge  of  angle 
2a  which,  to  simplify  the  presentation,  is  assumed  in  the  main  text  to  be 
such  that  7r/a  is  an  integer  m  (Figure  1).  The  wave  speeds  in  the  interior 
and  exterior  of  the  wedge  are  distinct  and  the  radiation  condition  of  only 
outgoing  waves  at  infinity  is  applied  in  all  directions.  At  the  boundary  of 
the  wedge  there  is  a  pair  of  transmission  conditions  which  ensure  continuity 
of  the  acoustic  pressure  and  normal  velocity.  Additional  simplification  is 
achieved  by  assuming  that  only  one  side  of  the  wedge  is  directly  ‘lit’  by 
the  source.  Such  a  transmission  problem  can  be  formulated,  by  use  of 
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the  free  space  Green’s  function  in  each  region,  in  terms  of  a  pair  of  coupled 
integral  equations  for  two  unknowns,  the  pressure  and  normal  velocity,  over 
the  total  boundary  of  the  wedge.  The  details  are  given  by  Kleinman  and 
Martin  (1988)  for  the  finite  transmitting  body. 

However,  by  using  suitably  modified  Green’s  functions  and  considering 
separately  the  symmetric  and  antisymmetric  parts  of  the  pressure  field  with 
respect  to  the  center  plane  of  the  wedge,  a  pair  of  disjoint  integral  equations 
of  the  first  kind  can  be  obtained  for  the  two  parts  of  the  normal  velocity 
on  just  one  face  of  the  wedge.  Transformation  to  equations  of  the  second 
kind  is  then  achieved  by  using  a  technique  for  solving  integral  equations 
with  Hankel  function  kernels  (Porter,  1983).  The  procedure  described  here 
can  be  regarded  as  a  perturbation  from  the  case  of  the  impenetrable  hard 
wedge  and  this  context  is  adopted  for  the  sake  of  extra  simplification  of  the 
presentation. 


3 


2  Formulation  of  the  Transmission  Problem 


A  sound  wave  generated  by  a  z-directed  source  of  period  2^/ uj  is  incident  on 
a  penetrable  wedge  which,  in  terms  of  cylindrical  polar  coordinates  (r,  <f>,  z), 
occupies  the  region  r  >  0,  —a  <  (f>  <  a  (Figure  1).  With  the  time  factor 
e'"‘  suppressed,  the  induced  exterior  and  interior  pressure  fields,  «e  and  «,• 
respectively,  satisfy  the  wave  equations 


{V^  +  kl)u,  =  0, 

(1) 

(V^  +  kj)u{  =  0, 

(l*^!  <  «). 

(2) 

while  the  incident  field,  with  the  source  at  {R,  $),  is  given  by 

Uinc  =  -^R?  - 

2rRcos((f>  - 

(3) 

where  a<$<7r  —  a,  in  accordance  with  the  assumption  that  only  one 
side  of  the  wedge  is  directly  ‘lit’  by  the  source.  The  transmission  conditions 
at  the  interfaces  are 


where 


u  =  Ui, 


du  du{ 


{<t>  =  ia)) 


U  =  +  Wine 


(4) 

(5) 


denotes  the  total  pressure  field  in  the  exterior  region.  At  infinity,  all  fields 
contain  only'  outgoing  waves.  The  wavenumbers  k{  and  the  density 
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ratio  p  are  given  real,  positive  constants  with  ki  >  k^.  Now  introduce  the 
symmetric  and  antisymmetric  parts  of  «e)  by  writing,  for  each  field 

u^^\r,(l>)  =  ^[u{r,<l))  ±u{r,-<f>)]  (6) 


and,  similarly,  set 

(>*>0).  (7) 

Then  the  exterior  and  interior  regions  are  reduced  to  a  <  ^  <  tt  and 
0  <  ^  <  a  respectively,  with  the  differential  equations  (1)  and  (2)  satisfied 
therein  and  the  transmission  conditions  (4)  applicable  at  the  sole  interface 
<j}  =  a.  The  additional  hard  and  soft  conditions  are 


IT 


=  0, 


^  =  0  {<i>  =  ’t), 


^  =  0,  «S->  =  0  (^  =  0).  (8) 

Let  G<*>(r,^;r',.^')  (a  <  <  tt)  be  two  Green’s  functions  defined 

in  the  exterior  region  where  they  satisfy  (1)  except  at  the  singularity  in 
the  prescribed  term  ^iH^\kc{r^  +  r'^  —  2rr'cos(^  —  Similarly,  let 

<  O')  be  two  Green’s  functions  defined  in  the  inte¬ 
rior  region  where  they  satisfy  (2)  except  at  the  singularity  in  the  prescribed 
term  ^iH^\ki{r^  +  r'^  —  2rr'cos(^  —  These  four  functions  satisfy 

the  respective  hard  and  soft  conditions  (8)  and  all  of  them  are  required  to 
have  zero  normal  derivative  at  the  interface  <j>  =  a  and  only  outgoing  waves 


at  infinity. 
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Then  when  Green’s  theorem  is  applied  to  <f>)  —  <j>)]p=Q  and 

r',  a)  in  the  external  region  a  <  <j>  <  ir  and  to  d)  and 

gS*V.  r',  a)  in  the  interior  region  0  <  ^  <  a,  the  resulting  pairs  of 
integral  equations  on  the  wedge  boundary  are 


«)  -  “)]p=o  =  -  /  a;  /,  a)  dr' 

TT  JQ  T  0(p 

=  —  f  a;  J*',  a)f^'^\r')dr'  (r  >  0) 

TT  Jo 

and 

u\^\r,  a)  =  G^^\r,  a;  /,  a)  dr' 

=  — —  /  a;  r',  a)/^'^^(r')dr'  (r  >  0), 

after  use  of  (4),  (5)  and  (7).  Now  (4)  and  (5)  allow  the  unknown  functions 
on  the  left  hand  sides  to  be  eliminated  to  obtain  a  pair  of  disjoint  integral 
equations  of  the  first  kind  for  the  symmetric  and  antisymmetric  source 
density  functions  namely 


«)  +  “)]p=o  =  a)]p=o 

=  — ^  j  a;  r',  a)  +  pG^^\r,  a;  r',  a)]  f^^\r')dr'  (r  >  0). 

(9) 

This  is  the  equation  whose  reduction  to  a  suitable  integral  equation  of  the 
second  kind  is  the  purpose  of  this  work.  Note  that  />  =  0  for  the  hard  wedge, 
in  which  limit  /  becomes  irrelevant  to  the  scattering  problem  because  no 
transmission  occurs  in  this  case.  However,  the  structure  of  the  solution 
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procedure  presented  below  may  be  regarded  as  a  perturbation  about  this 
limit.  On  writing 

=  +  +  + .  (lO) 


the  integral  equation  (9)  reduces  to  the  system 


a)]p=o  =  “”  / 

TT  •/O 

i  r  G<*>(r, «;  r',  a)fi^^(r')dr'  =  --  /”  gS^V.  “)/?’(■•')*'  (n  >  1) 

IT  Jo  ’T  JO 


(r  >  0). 


(11) 


Thus  each  term  in  the  expansion  (10)  is  determined  from  an  interior  prob¬ 
lem  while  each  forcing  term,  after  the  first,  arises  from  an  exterior  problem. 
In  particular,  this  shows  that  the  Green’s  function  with  wavenumber 
ki  plays  the  dominant  role  in  the  kernel  of  (9). 


3  The  Green’s  Functions  and  Forcing  Term 


With  the  incident  field  defined  by  (3),  it  may  be  deduced  from  Jones  (1986, 
section  9.19)  that  the  forcing  term  in  (9)  is  given  by 


4±)(r)  =  [u^^\r,a)],=o  =  - T  HS‘\k^{r^+R^-2rR  cosh 

'  lO^TT  —  aj  Joo-7rt 

1 


X  sinh 


ITT 


2{v  —  a) 


L<=osh5^-cosf5g  coshj^-hcosgiig 


dr 


=  +  e-  2rJ!cos(*  -  a))''^] 
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+— - r  /  [ke{r^  +  R^  +  2rR  cosh  u)^/^]  sinh 

lo(7r  —  a)  J-oo  /(tt  —  a) 

cosh  —  cos  ^  cosh  +  cos 

in  which  only  the  first  term  contributes  to  the  residue  at  r  =  0  because  of 
the  assumption  $  +  a  <  tt. 

Similarly  the  external  Green’s  functions  are  given  by 

1  roo+fft  HQ^\ks{r^  +  r'^  —  2rr'  cosh  sinh 


,  V  1  /OO  +  TT 

Gi+)(r.  <i>\  r',  a)  =  ^  /  . 

O^TT  —  O' j  ^oo— TTi 


cosh  +  cos 


dr, 


i.e. 


Gi+V,  «;<•',  a)  = 


+ 


8( 


— r 

TT  —  Of)  d-c 


“  +  2rr'cosh 


tanh#±2«l 
2(7r— a) 


du 


(13) 


and 


,  V  Sin  Tf - -ir  foo- 

4orr^r/=- 


=+ir.'  +  r'^  -  2rr'coshr)^/*]  sinh 


cosh  +  cos 


dr 


i.e. 


Gl->(r,a;r',a)  =  -i^r[Ae|r-r'|] 


+ 


_J_r 

(tt  —  a)  J-c 


°°  Ho^\ke{r^  +  r'^  +  2r/ cosh 


sinh^^ 

2(7r— a) 


du. 


(14) 


8(7r  —  a) 

The  assumption  that  a  =  for  some  integer  m  >  1,  simplifies  the 

corresponding  expressions  for  the  interior  Green’s  functions  which  evidently 
can  be  constructed  by  adding  image  sources  and  sinks  to  the  prescribed 
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source.  Thus 


-1  n=m 

<j)-,  r',  O')  =  -i  ^  (±1)"  +  r'^  —  2rr'cos[(2n  -  1)q;  + 

+  r'^  —  2rr'cos[(2n  —  1)0;  — 

and,  in  particular, 

«;  »•'>  “)  =  -  *•'1]  +  (±l)"*^o^^[*,(r  +  r')] 

n=m—l  ^ 

+2  ^  (±i\Y H^\ki{r^  +  r'^ —  2rr' cos .  (15) 

n=l  J 

4  Inversion  of  the  Hankel  Kernel 


Consider  the  integral  equation  of  the  first  kind 


g[r)  =  /  il){r')Hf^{ki\r  -  r'\)dr'  (r  >  0),  (16) 

^0 


whose  solution,  which  may  have  an  integrable  singularity  at  x  =  0,  is 
assumed  to  be  bounded  at  infinity,  and  define  the  operator  X,-  by 


{Lii)){s)  =  -ij)'{s)  +  kil  ^(r) 

Jg  —  5 

Porter  (1983)  showed  that  application  of  this  operator  to  (16)  yields  the 
singular  integral  equation 


Oi  too  pilti(T'-s) 

{U9){s)  =  -f  /  {s  >  0), 

TT  JO  r'  —  5 


(18) 


9 


whose  solution  is  given  by  Muskhelishvili  (1953).  The  complementary 
function,  /\/r',  satisfies  the  outgoing  wave  condition  and  furnishes 

a  unique  solution  bounded  at  the  origin,  namely 

m  r  (r  >  0).  (19) 

2t  Jo  y/s  s  -  r 

Thus,  on  writing  (15)  in  the  form 

gS*V.  >•',  a)  =  i.-  -  --'I)  +  JfPV.  >■')]  ,  (20) 

the  first  of  the  integral  equations  (11)  can  be  rearranged  as 

4Tia<'‘>(r)  -  r  ’•')/$*’(>■')*'  =  r 

Jo  Jo 

which  is  of  the  type  (16)  when  the  left  hand  side  is  regarded  as  known. 
Hence,  from  (19), 

=  Vrj:  [2(£,.S^-)(.)  +  (£<l^‘')(./)#>(/)^/] 

(r  >  0),  (21) 

which  is  an  integral  equation  of  the  second  kind  for  each  zero  order  density 
function.  The  two  integrations  introduced  by  the  above  inversion  procedure 
can  be  evaluated  exactly  when  plane  wave  representations  are  used. 

The  inversion  of  the  remaining  integral  equations  in  (11)  can  be  achieved 
similarly,  except  that  the  forcing  terms  also  have  a  Hankel  kernel,  with 
wavenumber  as  given  by  (13)  and  (14)  for  the  symmetric  and  antisym¬ 
metric  cases.  On  writing  these  equations  in  the  form 

G<«(r,  a;  r',  a)  =  ji  [jfi”(i.|r  -  r'|)  +  Ki*\r,  r')]  ,  (22) 


ds 
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as  in  (20)  for  the  interior  Green’s  functions,  the  sequence  of  integral  equa¬ 
tions  in  (11)  can  be  arranged  as 

-  (’•')  + 

-  r  w>(*.|r  -  r'l)  -  (fcl--  -  >-'l)]  fit\(r')dr' 

JQ 

=  HS‘\ki\r  -  r'l)  [/i^)(r')  +  fit\{r')]  dr', 
which  also  can  be  regarded  as  of  the  type  (16).  Hence,  on  defining  the 


bounded  kernel 


K,{r,  r')  =  -  r'|)  -  -  r'|)  (r,  r'  >  0),  (23) 

the  inversion  formula  (19)  yields 

^  f  f  /)£!(/) 

+{LiK.)is,  r')fit\{r')  -H  r')fi^\r')]  dr'ds  {n  >1)  (r  >  0), 

(24) 

which  is  an  integral  equation  of  the  second  kind  for  f^\n  >  1)  when  /^_i 
is  already  determined. 

The  expansion  (10)  enables  equations  (21)  and  (24)  to  be  combined  to 
show  that  the  corresponding  inversion  of  (9)  is 

(i+rt/<*)W  =  yr^  +  ^  /„ 

+piLiKi*''){s,T')  +  p(liK.)U,r')}f^’^\T')dT']ds  (r>0),  (25) 

which  is  an  integral  equation  of  the  second  kind  for 
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5  Evaluation  of  the  Inversion  Integrals 


The  two  integrations  introduced  by  the  above  inversion  procedure  can  be 
evaluated  exactly  when  plane  wave  representations  are  used.  Consider  the 
inversion  of  for  real  values  of  jS.  First,  it  is  readily  shown  from  the 


definition  (17)  that 


(i,e  )l,s)  -  I  isgn^(/92  - 


(1/91  <  *.  ) 
(I/9|  >  *.) 


Second,  the  evaluation  of 


o«fc,<5-r) 


^  Jo  y/s(s  —  r) 

requires  that  of  the  integral  /  given  by 


I(\,r)  =  Vrl 


00 


-ds  (r  >  0), 


-'Jo  Vs{s-rr  " 

which  vanishes  at  A  =  0  and,  by  consideration  of  can  be  expressed  in 
terms  of  the  Fresnel  integrals  Ci  and  Si  defined  by 


Ci{x)  +  iSiix)  =  ^-  r  e'^dy. 

T  Jo 


Thus 


J  =  -irv'2e-'''‘  {Cl  [(Ar)‘'“]  -  tSiKAr)*'"]}  (A  >  0), 

/  =  -,rv'2e''''*{C,(C|A|r)‘'^]  +  iS.[(|Alr)‘'=l}  (A  <  0).  (29) 

Hence  it  follows  from  (26)  and  (29)  that 

■Jr  r  r),  (30) 

Jo  ^  S  V  j 


12 


where 


Miiff,  r)  =  I  -(i?  -  /J2)‘/2e-W4  (CiCK*,-  -  «>•]"")  -  iSidCi.-  -  ^)rl‘/=)} 

(31) 

Now  in  order  to  use  (30),  it  is  necessary  to  use  plane  wave  representa¬ 
tions  of  the  functions  to  be  inverted  on  the  right  hand  side  of  (25).  The 
standard  representations  of  Jq  and  Iq  (Abramovitz  and  Stegun,  1964)  are 
combined  in  the  integral  representation 

H^^\k\r  -  r'l)  =  -  f  (32) 

TT  Jc 

where  C  is  the  contour  from  — ioo  to  tt  -f  icxi  drawn  along  the  negative 
imaginary  axis,  the  real  axis  from  0  to  tt  and  a  line  parallel  to  the  imaginary 
axis.  Equation  (32),  which  allows  (18)  to  be  derived  from  (16),  can  be 
generalized  to 

H^^\k{{x  -  x'f  +  {y-  (33) 

TT  Jc 

Similarly,  the  integral  representation 

ir?^(ifc|r  -  r'l)  =  r 
TTZ  J—oo 

can  be  generalized  to 

Hi^\k{{x  -  x'Y  -{y-  y'YY^^]  =  — ^ 

TTZ  */-oo 

(N  -  *1  >  »  -  »')•  (34) 


(^  >  *i) 

(||9|  <  *i) 
iff  <  -ki) 
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Now  the  forcing  and  kernel  functions  in  (25)  can  be  reduced,  as  follows, 
to  single  or  rapidly  convergent  double  integrals  involving  at  most  the  Fresnel 
integrals.  The  function  Uo*^(r),  given  by  (12),  can  be  written  in  the  form 


U'o 


(*)(r)  =  +  R^-  2rR  cos($  -  a))'/"*] 


4 

+  \  f  HQ^\ke{r^  +  R^ +  2rRcoshuY^^]h^'^\u)du 

2i  J  —  OQ 

_  _i_  j  g-»fc.[rcoBr-/icos(#— a+r]^^ 

47r  Jc 

k<'^\u)du  r  e-*‘[rco.h»+fico.h(«+v)]^y^  (35) 

27ri  J—oo  J—oo 


by  use  of  (33)  and  (34),  and  hence,  from  (30), 

2y/r  r  f^"'\{LiXif'>)(s)ds  =  -i-  / 

*/o  “*•  Tj  ^ 

+  iV2  r  h^^\u)du  r  Mi{K  cosh  u,  r)e-‘*‘t'  «+fico.h(«+v)]^y  (3g) 

J—OO  J—oo 

where  Mi  is  defined  by  (31).  The  contribution  from  the  interior  Green’s 
function  arises,  according  to  (15)  and  (20),  from 

n=m— 1 

Kj^\r,r')  =  (±l)”‘4"^[A.(r+rO]+2  £  {±1)’^ HS‘\ki{r^+r'^-2rr' cos  nz/myf^] 


— Tfl. —  J. 

g-.fc^rcosr  ^  (±l)ng-.V'co.(n:r/m+r)^^^ 

'  11=1 

by  use  of  (33).  The  first  term  can  be  inverted  by  comparison  with  (18)  and 
hence,  after  further  use  of  (30), 


ly/r  /■“ 

25r  Jo 


gifci(<-r) 
v/s(s  —  r) 


(LiK^*^){s,r')ds 


(±1)”* 

TT  \rV  r  +  r' 
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.  n=m— 1 

+  iy/2  COST,  r)  Y1  (37) 

n=l 

The  contribution  from  the  exterior  Green’s  function  arises  from  given 
by  (13),  (14)  and  (22),  and  K„  defined  by  (23).  On  using  (34)  and  then 
(30),  it  follows  that 

i^r  f°° 


tyr  /■= 

27r  Jo 


■\/s{s  —  r) 


r')*  = 


2v^27r(7r  —  a) 


/OO  Hfti  fOO 

- 7l/.(ifceCoshu,r)e-‘'"t'-'°=»'''+'''=°'>^(“+’')ldz;.  (38) 

‘OO  tanh  ■  ■  —  V  •/“OO 


2(7r— a) 


Similarly 


i^r 


I 


27r  Jo  y/s{s  —  r)  ‘  ®  ’  2\/27r(7r  —  a) 

r  ~ — ^J(u+Vi)  r  M.(i!:eCosht;,r)e-‘**f'~“‘’’'+''‘’“‘'("+’')ldu.  (39) 

J-oo  sinh  y ■■  -X  •'-°o 

ZylT—Ct ) 

For  the  remaining  term  in  (25),  write 


(LiK.){s,  rO  =  -  r'i))(»,  r')  -  -  r'|)](,,  r') 

-U,  -  U)H'^\K\r  -  r'|)](s,  r') 

2i 


It 


r'  —  s 


-  [(£.  -  Li)Hf(k.\r  -  r'|)]{»,  r')  (»,  r'  >  0), 

as  in  the  derivation  of  (18).  The  contribution  of  the  first  term  to 
iy/r  t°° 


(40) 


r 

27r  Jo 


y/s{s  -  r) 


iLiK.){s,r')ds 
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is  therefore 


Jo 


TT^ 


v/s(s  —  r)(r'  -  s) 


ds 


1  QtHr  -r)  /  r\^l^ 

=  - - [/(*.•  -  K, >•)-(-)  nk  -  k., r')l. 

r'  —  r  \r J 

where  I  is  defined  by  (27).  For  the  second  term  in  (40),  rewrite  (32)  as 


-  r'|)l(.,  r')  =  i  /_' 

vJi  (v^  -I 


'(t;2_  1)1/2 

and  observe  that  the  formula  for  {L^  —  Li)  corresponding  to  (26)  is 


(41) 


pe-i:.)e-‘''"”1(s)  =  e-'"'- 


ke{l  —  —  {kf  —  klv^y/^  (|t;|  <  1) 

isgnvk,{v^  -  1)1/2  _  (Jt?  _  ^2„2)l/2  (1  <  |„|  < 

isgnu[A:e(u2  —  1)^/2  —  {klv^  —  A:2)i/2]  (|u|  >  ki/ki) 

(42) 


The  factors  on  the  right  hand  side  of  (42)  effectively  replace  those  in  (26) 
and  account  must  be  taken  of  this  when  adapting  (30)  in  order  to  use  the 
function  M{  given  by  (31).  Hence  the  above  results  imply  that 

1  /  r  \  1/2 


roo  pilfi(M-r)  1  gt/iii’’  -r)  .  /  T  \  1/2  ,  ' 

/  47 - AL,K.){,,r')ds  =  -,—, - /ti-t.,/ 

Jo  ^ys{s-r)^  TT^  r'-r  '  \rJ 

+-^  r  e-'"('-'>M.p,u,r) 

TTy/l  J-1 


[{kf  —  kiv^yi^  (1  — 

-  rie-'X^-^'^Miik.v,  r)tp{v)  +  r)i^{-v)]dv,  (43) 

TT  '^r  2  J 1 

where  I  is  given  by  (29)  and 

(fc?-fcii;2)'l/2  -  („2_\)l/2  (1  <  |ll|  <  h/ke) 

~  (v2- 1)1/2  (ll/l  >  ^i/K) 


lj){v) 


(44) 
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Evidently  the  expression  (43)  is  bounded  at  r  =  r\  with  the  integral  re¬ 
maining  convergent  in  this  limit  because  the  integrand  is  0{v~^)  as  u  — »•  oo. 

Thus  equations  (36-39)  and  (43)  furnish  simplified  forms  of  the  integrals 
on  the  right  side  of  the  integral  equation  (25).  Additional  similar  terms  will 
be  introduced  on  relaxation  of  either  or  both  of  the  assumptions  a  =  7r/2m 
and  $  -h  a  <  ir.  Asymptotic  estimates  for  large  r  show  that  must  be 
a  linear  combination  of  and  ^  limit,  as  might  be 

anticipated.  This  information  will  facilitate  a  numerical  solution  for  the 
even  and  odd  components  of  the  normal  velocities  at  the  interfaces  of  the 
wedge  of  angle  2a. 
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Abstract.  A  pair  of  coupled  integral  equations  of  the  second  kind  are  specialized  to  the  even 
and  odd  symmetry  components  for  the  scattering  by  the  penetrable  wedge.  The  unknowns 
are  both  the  scalar  field  emd  its  normal  derivative  on  the  two  wedge  surfaces  that  separate 
the  interior  and  exterior  regions  of  different  wave  speeds  and  constitutive  parameters  (den¬ 
sity  or  dielectric  constant,  depending  on  the  specific  physical  application).  A  discretization 
and  collocation  procedure  transforms  the  operator  equations  into  coupled  matrix  equations, 
which  are  solved  numerically.  Further  extensions  and  plans  to  include  anticipated  asymptotic 
behavior  far  from  the  line  source  and  wedge  apex  are  outlined. 


I.  Adaptation  of  Kleinman  and  Martin  Analysis 

The  transmission  problem  for  the  penetrable  wedge  of  Figure  1  consists  of  appropriate 
exterior  and  interior  waves 

+  kl)  Uf.{r,<f>)  = -^8{r  —  ro)6(<j)  -  <f>Q)  {a<<f><2ir-a)  (1) 

(V^  +  kj)  Ui(r,  ^)  =  0  {-a  <  (f>  <  a)  (2) 

subject  to  the  two-dimensional  Sommerfeld  radiation  condition  plus  the  continuity  bound¬ 
ary  conditions 

Ue{r,±a)  =  Ui{r,±a)  (3) 

d  0 

— Ue(»*,  ia)  =  p—Ui{r,  ia)  (4) 

where  p  =  pel  Pi  is  the  ratio  of  exterior  to  interior  constitutive  parameters  (ambient  density 
or  dielectric  constant).  A  pair  of  integral  equations  derived  by  Kleinman  and  Martin  [4] 
for  the  scalar  field  u  and  its  normal  derivative  du/dn  on  the  boundary  of  a  penetrable 
obstacle,  are  adapted  to  the  semi-infinite  wedge  faces  of  Figure  1  to  effect  an  accurate 
numerical  solution  of  the  scattering  problem. 


Typeset  by 


I 


Fig.1.  Penetrable  Wedge  and  Line  Source 


Define  the  surface  distributions 


v±{r)  =  «(r,  O')  ±  «(r,  —a) 


«^±(r)  =  ■^u{r,a)  ±  -a) 


(5) 


(6) 


in  order  to  succinctly  decompose  the  field  into  its  symmetric  and  antisymmetric  compo¬ 
nents  with  respect  to  the  wedge  bisector  (the  x  axis). 

After  careful  algebra,  the  selected  integral  equations  from  [4]  are  written  in  a  form  de¬ 
fined  on  a  single  semi-infinite  domain  (0<r<oo)to  account  for  the  even/odd  symmetry: 


dr'  [G'e(r,Q';r',-a)  -/5G',(r,Qr;r', -O')] 

0 

OO 

-  J  dr'  u;±(r'){[Ge(r,  a;  r' ,  a)  -  G,(r,  a;  r' ,  a)]  ±  [G«(r,  a;  r',  -a)  -  G,(r,  a;  r',  -or)]} 

0 

=  2v^%r)  (7) 


OO 

/Q 

dr'  w±{r')—  [pGe{r,a-,  r',-a)  -  G,(r,  a;  r',-a)]  +  p 


■  J  dr'  v±{r') 


^2 

''  { [Geir,  a;  r',  a)  -  G.  (r,  a;  r' ,  a)]  ±  [Ge{r,  a;  r',  -a)  -  G.(r,  a;  r',  -a)] } 

=  2pu;r(r).  (8) 


2 


This  arrangement  ensures  regular  or  (at  worst)  weakly  singular  kernels,  which  involve 
the  dilFerences  of  the  scaled  interior  and  exterior  free-space  Green’s  functions  and  normal 
derivatives  thereof: 


[G.(r,  a;  /,  -a)  -  pG,{r,  a;  r',  -a)\  =  [pk,H[»{kR)  -  (9) 


dn 


Geir,a-,r',a)-Gi{r,a]r',a)  =  —  H^^\ke\r  -  r'\)  -  H^^\ki\r  -  r'\) 


-\n^ 


t-*t‘  TT  ki 

G,(t,  a;  r',  -a)  -  Gi(r.  a;  r',  -a)  =  i  J*)] 


dn 


2i 

[pG.ir,  a;  r',  -a)  -  Gi(r,  a;  r',  -a)]  =  -  pk.n['Hk,R)\  (12) 


(10) 


(11) 


dndn 


-  [Ge{r,  a;  r',  or)  -  Gi(r,  a;  r',  a)] 


-  r'l)  -  kiHP{ki\r-r'\)] 


2i\r  —  r'l 

ke\r-r'\ 


- . '-i-fe?ln 

r— ►r'  27r  \ 


-) 

/(13) 


dndn 


-  [Ge{r,  a;  r',  -a)  -  G,-(r,  a;  r',  -a)] 


2iR  {  ' 


cos  2a  — 


2rr'  sin^  2a 


kiH['\k,R)  -  k,Hl'\k,R)\ 


[i,?fr<'>(l;fll)  -  klH^'\k.R)\  I .  (14) 


The  distance 

R={r^+r'^-2rr'cos2a)^f^  (15) 

between  points  on  opposing  wedge  faces  is  zero  only  as  both  points  coalesce  at  the  apex 
(r  =  r'  =  0).  Excitation  from  the  line  source  of  unit  strength  at  (ro,  ^o)  gives  the  forcing 
terms 
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ikeVQ 


sm{(f>o  -  O') 


(jCeV'"'^  +  j’o  “  2rrocos(^o  —  o:)^ 

\/r^  +r^  -  2rro  cos(^o  -  «) 

h[^'^  (jces/r^  +  j’o  “  2rro  cos(^o  +  «)) 


±  sin(«^o  +  a) 


v/r^  +  Tq  —  2rro  cos{<f>o  +  a) 


(17) 


II.  Piecewise  Constant  Approximation  with  Collocation 
A  moment-method  expansion  in  terms  of  the  pulse  functions  of  width  A  in  r, 

1,  (m  —  1)A  <  r  <  mA 
0,  otherwise 

proceeds  formally,  without  regard  to  truncation  or  anticipation  of  asymptotic  behavior,  as 


v±(r)  =  ^  v^Pm(r),  w^Pm{r).  (19) 

iCp 

m=:l  m=l 

Collocation  at  the  center  points  of  each  pulse 

r,  =  (/-l/2)A  (20) 

yields  the  coupled  matrix  equations 

(1  +  ±  MI-'*]  -  [B*!!®*]  =  [/*)  (21) 

(1  +  rt[/l[«;±]  T  [Cl[u.*l  +  =  p\g^].  (22) 

Upon  discretization,  the  attractive  feature  of  integral  equations  of  the  second  kind  is  main¬ 
tained  in  the  form  of  diagonal  dominance.  The  elements  of  the  column  vectors  [/*]  and 
[g^]  are  the  right  hand  sides  of  (7)  and  (8),  respectively,  evaluated  at  the  collocation  points 

/f  =  24"'(r,),  gf  =  y/'M-  (23) 


and  [7]  is  the  identity  matrix. 

Integral  forms  for  the  dimensionless  coefficients  are 

mA 

Aini=  J  dr' -^[Ge{ri,a-,r',-a)  -  pGi{ri,a\r',-a)]  (24) 

(m  — 1)A 
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r 


mA 

B^=ke  j  dr' {[Geiri,a-,r',a) -Gi{rua-,r',a)] 

(m  — l)A 

±[Ge(r/,a;r',-Q;) -Gi(rj,a;r',-a)]}  (25) 

mA 

Cim  =  j  [pGe{ri,a]  r',  -a)  -  G,(r/,  a;  /,  -a)]  (26) 

(m  — 1)A 

mA 

(m  — 1)A 

±  [G«(ri,Q!;r',-a)  -  G,(rj,  a;  r', -a)]}.  (27) 

III.  Computational  Details 

Machine  computation  is  facilitated  by  introducing  the  normalized  parameters 
x' =  k^r'  xi  =  keri  (r  =  ki/ke 

8  =  keA  F  =  keR  =  {x]  +  x'^ -2xix' cos2aYf^.  (28) 

The  required  matrix  elements  are  now  obtained  by  quadrature,  with  careful  attention  to 
extract  the  removable  and  logarithmic  singularities  that  occur  in  two  species  of  the  diagonal 
terms  (/  =  m): 


■^Im  — 


xi  sin  2q! 
2i 


mS 

J  dx'  j 


(29) 


(m  — 1)6 


pi  _  -1- 


mb 


B 


(1) 

Im 


=  / 


H. 


-  x'\)  -  Hi^^{<r\x,  -  i'|)' 


(m  — 1)6 


2i 


(30) 


•  -J*  In  as  x*^xi 


mb 


bS  =  ^  Vf )] 


c,„  = 


sin  2a 
2i 


(m  — 1)6 
m6 


(31) 


(32) 


(m  — 1)6 
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'/m 


Im 


mS 


«-  / 


,  -  x'\)  -  (rH'l^>{(r\xi  -  x^|) 

2i\xi  —  x'\ 


(1), 


dx 


(m  — 1)6 


(33) 


■  ^In  ^*^2*  ^  as  c'— ►ajj 


i) 


(2) 

Im 


m6 

['“2a- 


2xix'  sir?  2oi 


F2 


affPVf)  -  ffi'V )] 


(m— 1)5 


+ 


xix'  sin^  2a 


iT^B^iirF)  -  )]  }.  (34) 


A  Fortran  program  that  implements  the  above  analysis  is  included  in  Appendix  A.  A 
sharp  truncation  of  both  the  expansions  (19)  and  the  collocation  points  (20)  a,t  l,m  =  N 
results  in  close  (not  large  r)  behavior  of  t;(r)  and  w{r)  that  is  surprisingly  persistent. 
However,  an  accurate  solution  must  account  for  the  infinite  domain.  A  physically-based 
approach  is  to  include  both  cylindrical  waves  that  Davis^  extracts  from  an  asymptotic 
estimate: 

N 


w±{r) 


N 


g.fcir  ^  g.fc.r 

(35) 

^ikiT  ^  ^ik^T 

(36) 

The  infinite-range,  highly  oscillatory  integrals  required  by  this  modification  are  calculated 
with  the  aid  of  a  x/2  rotation  in  the  complex  plane,  resulting  in  exponentially  decaying 
integrands  that  are  more  suitable  for  direct  quadrature. 


IV.  Extensions  and  Conclusions 

Accurate  and  efficient  implementation  of  the  above  mathematics  requires  careful  nu¬ 
merical  analysis.  The  proposed  asymptotic  ansatz  to  account  for  the  far  behavior  of  the 
unknown  surface  distributions  is  best  verified  first  for  a  simpler,  known  problem:  the  soft 
(or  hard)  half-plane.  Although  the  selected  set  of  coupled  integral  equations  does  intro¬ 
duce  a  pair  of  unknown  functions,  as  opposed  to  some  methods  that  require  only  a  single 
unknown,  the  attractive  features  are  weakly  singular  kernels  and  second-kind  structure. 
Also,  having  direct  access  to  both  the  scalar  field  and  its  normal  derivative  on  the  boundary 
is  a  benefit  to  the  techniques  of  far-field  evaluation  based  on  various  Green’s  functions. 


^A.M.J.  Davis,  Acoustical  Scattering  by  a  Penetrable  Wedge,  1994,  p.  17,  in  this  report. 
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Appendix  A.  Fortran  Program  for  Coupled  Integral  Equations 


o  o  on 


C  FILE  "COUPSEQ  FORTRAN  A"  -  COUPLED  INTEGRAL  EQUATIONS  FOR 
C  DIELECTRIC  WEDGE  SCATTERING. 

C  15  MAY  1994  R . W . SCHARSTEIN  UNIVERSITY  OF  ALABAMA  EE  DEPT. 

C 

PARAMETER  (NN=50) 

EXTERNAL  FA , FAREAL , FAIMAG , FB 1 , FB IREAL , FB 1 IMAG , FB2 , FB2REAL , FB2 IMAG 

EXTERNAL  FC , FCREAL , FCIMAG , FD2 , FD2REAL , FD2IMAG 

EXTERNAL  FD 1 A , FD lAREA , FD lA IMA , FD IB , FD IBREA , FD IB IMA 

COMMON  /COM/XL, ALPHA, RHO, SIGMA 

COMMON  /SOURCE/XO,PHIO,ALPHA1 

COMPLEX  A(NN,NN) ,B1(NN,NN) ,B2(NN,NN) ,C(NN,NN) ,D1(NN,NN) ,D2(NN,NN) 
COMPLEX  B(NN,NnLd(NN,NN) 

COMPLEX  Q(2*NN,2*NN) ,Y(2*NN) ,VINC,WINC,V,W 
COMPLEX  CR0UT(2*NN,4*NN) , TEMP, DET, BIG, ERROR 
DIMENSION  LP(2*NN) 

INTEGER  P 
C 

DATA  PI/3.1415927/ 

CHOOSE  EVEN  (P=2)  OR  ODD  (P=l)  SYMMETRY 
P=2 


COUOOOlO 

COU00020 

COU00030 

COU00040 

COU00050 

COU00060 

COU00070 

COU00080 

COU00090 

COUOOlOO 

COUOOllO 

COU00120 

COU00130 

COU00140 

COU00150 

COU00160 

COU00170 

COU00180 

COU00190 

COU00200 

COU00210 

COU00220 


QUADRATURE  PARAMETERS 
EREL  =  0.01 
EABS  =  0. 

C  PHYSICAL  CONSTANTS 
ALPHA  =  PI/ 6. 

ALPHA 1  =  ALPHA 
RHO  =0.1 
SIGMA  =  SQRT(10.) 

C  SOURCE  COORDINATES 
XO  =  2. 

PHIO  =  PI/2. 

C  MORE  INTEGRATION  PARAMETERS 
XPMAX  =  25. 

DELTA  =  XPMAX/FLOAT(NN) 

C 

DO  10  L=1,NN 

XL  =  (FL0AT(L)-0.5)*DELTA 

DO  10  M=1,NN 

XI  =  FLOAT(M-l)*DELTA 

X2  =  FLOAT(M)*DELTA 

CALL  QDAG ( FAREAL , X 1 , X2 , EABS , EREL , 6 , AR , EREST) 

CALL  QDAG ( FAIMAG , X 1 , X2 , EABS , EREL , 6 , A I , EREST ) 

A(L,M)  =  CMPLX(AR,AI)*XL*SIN(2.*ALPHA)/CMPLX(0.,2.) 
C  WRITE(6,100)  L,M,A(L,M) 

CALL  QDAG ( FB IREAL , X 1 , X2 , EABS , EREL , 6 , B IR , ERE  ST ) 

CALL  QDAG ( FB 1 IMAG , X 1 , X2 , EABS , EREL , 6 , B 1 I , EREST ) 
B1(L,M)  =  CMPLX(B1R,B1I) 

C  WRITE(6,100)  L,M,B1(L,M) 

CALL  QDAG ( FB2REAL , X 1 , X2 , EABS , EREL , 6 , B2R , EREST) 

CALL  QDAG ( FB2IMAG , X 1 , X2 , EABS , EREL , 6 , B  2 1 , ERE  ST ) 
B2(L,M)  =  CMPLX(B2R,B2I)/CMPLX(0.,2.) 

C  WRITE(6,100)  L,M,B2(L,M) 

CALL  QDAG( FCREAL , X 1 , X2 , EABS , EREL , 6 , CR , EREST ) 

CALL  QDAG ( FC IMAG , X 1 , X2 , E ABS , EREL , 6 , C I , EREST ) 

C(L,M)  =  CMPLX(CR,CI)*SIN(2.*ALPHA)/CMPLX(0.,2.) 

C  WRITE(6,100)  L,M,C(L,M) 

IF(L.EQ.M)  GO  TO  20 

CALL  QDAG ( FD lARE A , XI, X2, EABS, EREL, 6, DIR, ERE ST) 


COU00230 

COU00240 

COU00250 

COU00260 

COU00270 

COU00280 

COU00290 

COU00300 

COUO'0310 

COU00320 

COU00330 

COU00340 

COU00350 

COU00360 

COU00370 

COU00380 

COU00390 

COU00400 

COU00410 

COU00420 

COU00430 

COU00440 

COU00450 

C0U00460 

COU00470 

COU00480 

COU00490 

COU00500 

COU00510 

COU00520 

COU00530 

COU00540 

COU00550 

COU00560 

COU00570 

COU00580 

COU00590 

COU00600 


o  o 


CALL  QDAG(FD1AIMA,X1,X2,EABS,EREL,6,D1I,EREST) 

D1(L,M)  =  CMPLX(D1R,D1I) 

GO  TO  30 
20  CONTINUE 

CALL  QD AG ( FD IBRE A ,X1,X2,EABS, EREL , 6 , D IR , ERE  ST ) 

CALL  QDAG (FD1BIMA,X1,X2,EABS, EREL , 6 , D 1 1 , ERE  ST ) 

D1(L,M)  =  CMPLX(D1R,D1I)  +  DELTA*( ( 1 . -SIGMA**2)*(ALOG(DELTA/4 . ) 
&  -1.)  -  SIGMA**2*ALOG(SIGMA))/6. 2831853 

30  CONTINUE 

C  WRITE(6,100)  L,M,D1(L,M) 

CALL  QDAG ( FD2REAL , X 1 , X2 , E AB  S , EREL , 6 , D2R , ERE  ST ) 

CALL  QDAG(FD2IMAG,X1,X2,EABS,EREL,6,D2I,EREST) 

D2(L,M)  =  CMPLX(D2R,D2I)/CMPLX(0. ,2.) 

C  WRITE(6,100)  L,M,D2(L,M) 

10  CONTINUE 

BOTH  EVEN  (+)  AND  ODD  (-)  SYMMETRY  PROBLEMS 
DO  40  L=1,NN 
DO  40  M=1,NN 

B(L,M)  =  B1(L,M)  +  (-1)**P*B2(L,M) 

40  D(L,M)  =  D1(L,M)  +  (-1)**P*D2(L,M) 

C  THE  BIG  SUPER-MATRIX  Q 
DO  42  L=1,NN 
DO  42  M=1,NN 
Q(L,M)  =  (-1)**P*A(L,M) 

IF(L.EQ.M)  Q(L,M)  =  Q(L,M)  +  l.+RHO 
Q(L,M+NN)  =  -B(L,M) 

Q(L+NN,M)  =  RHO*D(L,M) 

Q(L+NN,M+NN)  =  (- 1)**(P+1)*C(L,M) 

IF(L.EQ.M)  Q(L+NN,M+NN)  =  Q(L+NN,M+NN)  +  l.+RHO 
42  CONTINUE 

C  THE  BIG  RIGHT-HAND  COLUMN  VECTOR 
DO  44  L=1,NN 
XL  =  (FLOAT(L)-0.5)*DELTA 
Y(L)  =  2.*VINC(XL,P) 

44  Y(L+NN)  =  2.*RHO*WINC(XL,P) 

C  INVERT  THE  MATRIX  Q  AND  SOLVE  THE  SYSTEM  OF  LINEAR  EQUATIONS 
C  USING  CROUTC  (INCLUDED) 

DO  46  L=1,2*NN 

DO  46  M=1,2*NN 

46  CROUT(L,M)  =  Q(L,M) 

CALL  CROUTC ( 2*NN , 2*NN , 0 , CROUT , 1 . E - 7 , DET ,  lERR ,  LP  ) 

WRITE (6, 700)  lERR 

700  FORMAT(//,7X, ’SUBROUTINE- "CROUTC"  lERR  =  ’,13) 

C  CHECK  MATRIX  INVERSION 
BIG  =  CMPLXCO. ,0.) 

DO  57  K=1,2*NN 

DO  57  L=1,2*NN 

TEMP  =  CMPLX(0. ,0  .  ) 

DO  58  M=1,2*NN 

58  TEMP  =  TEMP  +  CROUT(K,M+2*NN)*Q(M,L) 

ERROR  =  TEMP 

IF(K.EQ.L)  ERROR  =  TEMP  -  CMPLX(1.,0.) 

IF(CABS(ERROR) .GT.CABS(BIG))  BIG  =  ERROR 
57  CONTINUE 

WRITE(6,710)  BIG 

710  FORMAT(/,7X, 'BIGGEST  ERROR  IN  MATRIX  INVERSION  =  ' ,2(E12.5,2X)) 
C  RESULTANT  COLUMN  VECTORS 
WRITE(6,200) 

200  FORMAT(//,4X, 'N* ,14X, 'V' ,22X, 'W' ,/) 


COU00610 

COU00620 

COU00630 

COU00640 

COU00650 

COU00660 

COU00670 

COU00680 

COU00690 

COU00700 

COU00710 

COU00720 

COU00730 

COU00740 

COU00750 

COU00760 

COU00770 

COU00780 

COU00790 

COU00800 

COU00810 

COU00820 

COU00830 

COU00840 

COU00850 

COU00860 

COU00870 

COU00880 

COU00890 

COU00900 

COU00910 

COU00920 

COU00930 

COU00940 

C0U00950 

COU00960 

COU00970 

COU00980 

COU00990 

COUOIOOO 

COUOlOlO 

COU01020 

COU01030 

COU01040 

COU01050 

COU01060 

COU01070 

COU01080 

COU01090 

COUOllOO 

COUOlllO 

COU01120 

COU01130 

COU01140 

COU01150 

COU01160 

COU01170 

COU01180 

COU01190 

COU01200 


210  FORMAT(2X,I3,2(2X,E12.5,2X,F7.2)) 

DO  60  L=1,NN 

V  =  CMPLX(0.,0.) 

W  =  CMPLXCO. ,0.) 

DO  62  M=1,2*NN 

V  =  V  +  CROUT(L,M+2*NN)*Y(M) 

62  W  =  W  +  CROUT(L+NN,M+2*NN)*Y(M) 

WRITE(6,210)  L,CABS(V),CDEG(V),CABS(W),CDEG(W) 
60  CONTINUE 

100  FORMAT(2X,I3,2X,I3,2X,E14.7,2X,E14.7) 

STOP 

END 

C 

C 

c 

FUNCTION  FA(XP) 

COMMON  /COM/XL, ALPHA, RHO, SIGMA 
COMPLEX  FA,H1,H2 

F  =  SQRT(XL**2+XP**2-2.*XL*XP*COS(2.*ALPHA)) 

HI  =  CMPLX(BSJ1(SIGMA*F),BSY1(SIGMA*F)) 

H2  =  CMPLX(BSJ1(F),BSY1(F)) 

FA  =  (RH0*SIGMA*H1-H2)/F 
RETURN 
END 
C 
C 

FUNCTION  FAREAL(XP) 

COMPLEX  FA 

F AREAL  =  REAL(FA(XP)) 

RETURN 

END 

C 

C 

FUNCTION  FAIMAG(XP) 

COMPLEX  FA 

FAIMAG  =  AIMAG(FA(XP)) 

RETURN 

END 

C 

C 

c 

FUNCTION  FBl(XP) 

COMMON  /COM/  XL, ALPHA, RHO, SIGMA 
COMPLEX  FB1,H1,H2 
D  =  ABS(XL-XP) 

IF(D.LT. 0.001)  GO  TO  10 
HI  =  CMPLX(BSJ0(D),BSY0(D)) 

H2  =  CMPLX(BSJ0(SIGMA*D),BSY0(SIGMA*D)) 

FBI  =  (H1-H2) /CMPLXCO. ,2.) 

RETURN 

10  FBI  =  -ALOG(SIGMA)/3. 1415927 
RETURN 
END 
C 
C 

FUNCTION  FBIREAL(XP) 

COMPLEX  FBI 

FBIREAL  =  REAL(FB1(XP)) 

RETURN 

END 


COU01210 

COU01220 

COU01230 

COU01240 

COU01250 

COU01260 

COU01270 

COU01280 

COU0I290 

COU01300 

COU01310 

COU01320 

COU01330 

C0U01340 

C0U01350 

COU01360 

C0U01370 

COU01380 

COU0I390 

COU01400 

COU0I410 

C0U01420 

C0U01430 

COU01440 

COU01450 

COU01460 

COU01470 

COU01480 

C0U01490 

C0U01500 

C0U01510 

COU01520 

C0U01530 

C0U01540 

COU01550 

COU01560 

COU01570 

COU01580 

COU01590 

COU01600 

C0U01610 

COU01620 

C0U01630 

C0U01640 

COU01650 

COU01660 

COU01670 

COU01680 

COU01690 

COU01700 

C0U01710 

COU01720 

COU01730 

C0U01740 

COU01750 

COU01760 

COU01770 

COU01780 

C0U01790 

COU01800 


n  o 


c 

c 

FUNCTION  FBIIMAG(XP) 
COMPLEX  FBI 

FBIIMAG  =  AIMAG(FB1(XP)) 
RETURN 
END 
C 


FUNCTION  FB2(XP) 

COMMON  /COM/  XL, ALPHA, RHO, SIGMA 
COMPLEX  FB2,H1,H2 

F  =  SQRT(XL**2+XP**2-2.*XL*XP*COS(2.*ALPHA)) 
HI  =  CMPLX(BSJO(F),BSYO(F)) 

H2  =  CMPLX(BSJO(SIGMA*F),BSYO(SIGMA*F)) 

FB2  =  H1-H2 
RETURN 
END 
C 
C 

FUNCTION  FB2REAL(XP) 

COMPLEX  FB2 

FB2REAL  =  REAL(FB2(XP) ) 

RETURN 

END 

C 

C 

FUNCTION  FB2IMAG(XP) 

COMPLEX  FB2 

FB2IMAG  =  AIMAG(FB2(XP)) 

RETURN 

END 

C 

C 

C 

FUNCTION  FC(XP) 

COMMON  /COM/  XL, ALPHA, RHO, SIGMA 
COMPLEX  FC,H1,H2 

F  =  SQRT(XL**2+XP**2-2.*XL*XP*COS(2.*ALPHA)) 
HI  =  CMPLX(BSJ1(SIGMA*F),BSY1(SIGMA*F)) 

H2  =  CMPLX(BSJ1(F),BSY1(F)) 

FC  =  XP*(SIGMA*Hl-RHO*H2)/F 
RETURN 
END 
C 

c 

FUNCTION  FCREAL(XP) 

COMPLEX  FC 

FCREAL  =  REAL(FC(XP)) 

RETURN 

END 

C 

C 

FUNCTION  FCIMAG(XP) 

COMPLEX  FC 

FCIMAG  =  AIMAG(FC(XP)) 

RETURN 

END 

C 


COU01810 

COU01820 

COU01830 

C0U01840 

COU01850 

C0U01860 

COU01870 

COU01880 

C0U01890 

COU01900 

COU01910 

COU01920 

C0U01930 

COU01940 

COU01950 

COU01960 

COU01970 

COU01980 

COU01990 

COU02000 

COU02010 

COU02020 

COU02030 

COU02040 

COU02050 

COU02060 

COU02070 

COU02080 

COU02090 

COU02100 

COU02110 

COU02120 

COU02130 

COU02140 

COU02150 

COU02160 

COU02170 

COU02180 

C0U02190 

COU02200 

COU02210 

COU02220 

COU02230 

C0002240 

COU02250 

COU02260 

COUG2270 

COU02280 

COU02290 

COU02300 

COU02310 

COU02320 

COU02330 

COU02340 

COU02350 

COU02360 

COU02370 

COU02380 

COU02390 

COU02400 


c 

c 

FUNCTION  FDIA(XP) 

C  USE  THIS  ONE  FOR  OFF-DIAGONAL  (L  NOT=  M)  TERMS  (NO  SINGULARITY!) 
COMMON  /COM/  XL, ALPHA, RHO, SIGMA 
COMPLEX  FD1A,H1,H2 
D  =  ABS(XL-XP) 

HI  =  CMPLX(BSJ1(D) ,BSY1(D)) 

H2  =  CMPLX(BSJ1(SIGMA*D),BSY1(SIGMA*D)) 

FDIA  =  (H1-SIGMA*H2)/CMPLX(0. ,2.*D) 

RETURN 

END 

C 

C 

FUNCTION  FDIAREA(XP) 

COMPLEX  FDIA 

FDIAREA  =  REAL(FD1A(XP)) 

RETURN 

END 

C 

C 

FUNCTION  FDIAIMA(XP) 

COMPLEX  FDIA 

FDIAIMA  =  AIMAG(FD1A(XP)) 

RETURN 

END 

C 

C 

C 

FUNCTION  FDIB(XP) 

C  USE  THIS  ONE  FOR  DIAGONAL  (L  =  M)  TERMS  (WITH  SINGULARITY!) 

COMMON  /COM/  XL, ALPHA, RHO, SIGMA 
GOMPLEX  FD1B,H1,H2 
D  =  ABS(XL-XP) 

IF(D.LT. 0.001)  GO  TO  10 
HI  =  CMPLX(BSJ1(D),BSY1(D)) 

H2  =  CMPLX(BSJ1(SIGMA*D) ,BSY1(SIGMA*D)) 

FDIB  =  (H1-SIGMA*H2)/CMPLX(0. ,2.*D) 

FDIB  =  FDIB  -  (ALOG(D/2.)-SIGMA**2*ALOG(SIGMA*D/2.))/6. 2831853 
RETURN 

10  FDIB  =  CMPLX(0. ,0 . ) 

RETURN 

END 

C 

C 

FUNCTION  FDIBREA(XP) 

COMPLEX  FDIB 

FDIBREA  =  REAL(FD1B(XP)) 

RETURN 

END 

C 

C 

FUNCTION  FDIBIMA(XP) 

COMPLEX  FDIB 

FDIBIMA  =  AIMAG(FD1B(XP)) 

RETURN 

END 

C 

C 
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FUNCTION  FD2(XP) 

COU03010 

COMMON  /COM/XL, ALPHA, RHO, SIGMA 

COU03020 

COMPLEX  FD2,H1,H2,H3,H4 

COU03030 

F  =  SQRT(XL**2+XP**2-2.*XL*XP*C0S(2.*ALPHA)) 

COU03040 

HI  =  CMPLX(BSJ1(SIGMA*F),BSY1(SIGMA*F)) 

COU03050 

H2  =  CMPLX(BSJ1(F) ,BSY1(F)) 

COU03060 

H3  =  CMPLX(BSJO(SIGMA*F),BSYO(SIGMA*F)) 

COU03070 

H4  =  CMPLX(BSJO(F),BSYO(F)) 

COU03080 

Z  =  XL*XP*(SIN(2.*ALPHA))**2 

COU03090 

FD2  =  rcOS(2.*ALPHA)-2.*Z/F**2)*(SIGMA*Hl-H2)+Z*(SIGMA**2*H3-H4)/FCOU03100 

FD2  =  FD2/F 

COU03110 

RETURN 

COU03120 

END 

COU03130 

c 

C0U03140 

c 

COU03150 

FUNCTION  FD2REAL(XP) 

COU03160 

COMPLEX  FD2 

COU03170 

FD2REAL  =  REAL(FD2(XP) ) 

COU03180 

RETURN 

COU03190 

END 

COU03200 

c 

COU03210 

c 

C0U03220 

FUNCTION  FD2IMAG(XP) 

COU03230 

COMPLEX  FD2 

C0U03240 

FD2IMAG  =  AIMAG(FD2(XP)) 

COU03250 

RETURN 

COU03260 

END 

COU03270 

c 

COU03280 

c 

COU03290 

c 

COU03300 

FUNCTION  VINC(XP,P) 

COU03310 

COMMON  /SOURCE/  XO ,PHI0 , ALPHA 

COU03320 

COMPLEX  VINC,H1,H2 

COU03330 

INTEGER  P 

COU03340 

G1  =  SQRT(XP**2+XO**2-2.*XP*XO*COS(PHIO-ALPHA)) 

COU03350 

G2  =  SQRT(XP**2+XO**2-2.*XP*XO*COS(PHIO+ALPHA)) 

COU03360 

HI  =  CMPLX(BSJO(G1),BSYO(G1)) 

COU03370 

H2  =  CMPLX(BSJ0(G2),BSY0(G2)) 

COU03380 

VINC  =  (H1+(-1)**P*H2)*CMPLX(0. ,l.)/4. 

COU03390 

RETURN 

COU03400 

END 
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c 

COU03420 

c 

' 

COU03430 

c 

COU03440 

FUNCTION  WINC(XP,P) 

C0U03450 

COMMON  /SOURCE/  XO ,PHI0 .ALPHA 

COU03460 

COMPLEX  WINC,H1,H2 

COU03470 

INTEGER  P 

COU03480 

G1  =  SQRT(XP**2+XO**2-2.*XP*XO*COS(PHIO-ALPHA)) 

COU03490 

G2  =  SQRT(XP**2+XO**2-2.*XP*XO*COS(PHIO+ALPHA)) 

COU03500 

HI  =  CMPLX(BSJ1(G1),BSY1(G1)) 

COU03510 

H2  =  CMPLX(BSJ1(G2),BSY1(G2)) 

COU03520 

WINC  =  (SIN(PHI0-ALPHA)*Hl/Gl+(-l)**P*SIN(PHI0+ALPHA)*H2/G2) 

COU03530 

&  *CMPLX(0. ,X0)/4. 

C0U03540 

RETURN 

COU03550 

END 

C0U03560 

c 

COU03570 

c 

COU03580 

c 

C0U03590 

FUNCTION  CDEG(Z) 

COU03600 

c 

C  PHASE  ANGLE  IN  DEGREES  OF  A  COMPLEX  NUMBER  Z. 

C 

COMPLEX  Z 
ZI  =  AIMAG(Z) 

ZR  =  REAL(Z) 

IF((ZI.EQ.O.).AND.(ZR.EQ.O.))  GO  TO  10 
CDEG  =  ATAN2(ZI,ZR)*57. 29578 
RETURN 
10  CDEG  =  0. 

RETURN 

END 

C 

C 

C 

SUBROUTINE  CROUTC ( MR , NR , NCC , A , ZMGH , DT , lERR , LP ) 


C  GROUT  (1)  OPERATES  ON  A  COEFFICIENT  MATRIX  TO  SOLVE  A  SYSTEM  OF 
C  SIMULTANEOUS  EQUATIONS  OR  TO  COMPUTE  AN  INVERSE  AND  (2) 

C  COMPUTES  A  DETERMINANT. 

C  CROUT  REDUCES  THE  ORIGINAL  MATRIX  AND  RIGHT-HAND  SIDES  UNTIL 
C  UPON  COMPLETION  THE  REDUCED  MATRIX  REPLACES  THE  ORIGINAL 
C  MATRIX  AND  THE  SOLUTIONS  REPLACE  THE  RIGHT-HAND  SIDES. 

COMPLEX  *  8  A,DT,TEMPC 

DIMENSION  A(l),  LP(1) 

IF(NR.GT.MR)  GO  TO  210 
MTX  =  MR*NR 
MRA  =  MR  +  1 
MRS  =  MR  -  1 
MDN  =  MR  -  NR 
MTR  =  MTX  -  MDN 
MTRA  =  MTX  +  1 
DT  =  (l.,0.) 
lERR  =  0 
NRS  =  NR  -  1 
DO  2  I  =  1, NR 
LP(I)  =  I 
2  CONTINUE 

IF(NCC)210, 805, 1001 
805  NTC  =  NR  +  NR 

MTT  =  MR*NTC  -  MDN 
J  =  MTRA 

DO  19  K  =  MTRA, MTT, MR 
JF  =  K  +  NRS 
DO  18  KX  =  K,JF 
A(KX)  =  (0.,0.) 

18  CONTINUE 

A(J)  =  (l.,0.) 

J  =  J  +  MRA 

19  CONTINUE 
GO  TO  1 

1001  NTC  =  NR  +  NCC 

MTT  =  MR*NTC  -  MDN 
1  IF(NTC.LE.NR)  GO  TO  210 
DO  70  I  =  1,NR 
IS  =  I  -  1 
II  =  MR*IS  +  I 
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C0U04120 

COU04130 

COU04140 

C0U04150 

COU04160 

COU04170 

C0U04180 

COU04190 
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30 

33 


31 


34 

35 

45 

46 


48 

47 

70 

78 


160 

170 

180 


USB  =  11-1 

COU04210 

HAD  =  II  +  MR 

COU04220 

ICF  =  MR*I  -  MDN 

COU04230 

ICS  =  ICF  -  NRS 

C0U04240 

HA  =  11  +  1 

COU04250 

TEMP  =  0. 

COU04260 

DO  31  J  =  H,MTT,MR 

COU04270 

IF(I.EQ.l)  GO  TO  33 

COU04280 

KF  =  J  -  1 

COU04290 

KS  =  J  -  I  +  1 

COU04300 

KX  =  I 

COU04310 

DO  30  K  =  KS,KF 

COU04320 

A(J)  =  A(J)  -  A(KX)*  A(K) 

COU04330 

KX  =  KX  +  MR 

COU04340 

CONTINUE 

COU04350 

IF(J.GT.MTR)  GO  TO  31 

COU04360 

IF(CABS(A(J)).LE.TEMP)  GO  TO  31 

C0U04370 

TEMP  =  CABS(A(J)) 

COU04380 

NX  =  J/MR  +  1 

COU04390 

CONTINUE 

COU04400 

IF(I.EQ.NR)  GO  TO  35 

COU04410 

IF(NX.EQ.I)  GO  TO  35 

COU04420 

ITEMP  =  LP(NX) 

COU04430 

LP(NX)  =  LP(I) 

C0U04440 

LP(I)  =  ITEMP 

COU04450 

LPIS  =  MR*NX  -  MRS 

C0U04460 

DO  34  K  =  ICS, ICF 

COU04470 

TEMPO  =  A(K) 

COU04480 

A(K)  =  A(LPIS) 

COU04490 

A (LPIS)  =  TEMPO 

COU04500 

LPIS  =  LPIS  +  1 

COU04510 

CONTINUE 

COU04520 

DT  =  -  DT 

COU04530 

DT  =  DT  *  A(H) 

C0U04540 

IF(ZMCH  -  CABS(A(H)))  45,45,220 

COU04550 

DO  46  J  =  IIAD,MTT,MR 

COU04560 

A(J)  =  A(J)/A(H) 

C0U04570 

CONTINUE 

C0U04580 

IF(I.EQ.l)  GO  TO  70 

COU04590 

IF(I.EQ.NR)  GO  TO  78 

COU04600 

DO  47  M  =  HA, ICF 

COU04610 

KX  =  M  -  ICS  +  1 

COU04620 

DO  48  KY  =  ICS, USB 

COU04630 

A(M)  =  A(M)  -  A(KX)  *  A(KY) 

COU04640 

KX  =  KX  +  MR 

C0U04650 

CONTINUE 

COU04660  ' 

CONTINUE 

COU04670 

CONTINUE 

COU04680 

NRTAD  =  MTX  +  NR 

COU04690 

DO  180  I  =  1,NRS 

COU04700 

IREV  =  NRTAD  -  I 

COU04710 

KRS  =  IREV  -  MR* I 

COU04720 

DO  170  IRCNT  =  IREV,MTT,MR 

COU04730 

KCS  =  IRCNT  +  1 

C0U04740 

DO  160  K  =  KRS,MTR,MR 

COU04750 

A(IRCNT)  =  A(IRCNT)  -  A(KCS)  *  A(K) 

COU04760 

KCS  =  KCS  +  1 

COU04770 

CONTINUE 

COU04780 

CONTINUE 

COU04790 

CONTINUE 

COU04800 

DO  6  I  =  1,NRS 
9  IF(LP(I).EQ.I)  GO  TO  6 
NX  =  LP(I) 

LP  (I)  =  LP(NX) 

LP(NX)  =  NX 

IXS  =  MTX  +  I 

lY  =  MTX  +  NX 

DO  7  IX  =  IXS,MTT,MR 

TEMPO  =  A(IX) 

A(IX)  =  A(IY) 

A(IY)  =  TEMPO 
lY  =  lY  +  MR 
7  CONTINUE 
GO  TO  9 
6  CONTINUE 
RETURN 

210  lERR  =  2 
NTC  =  NR 

DT  =  (999999999 . ,999999999 . ) 
RETURN 

220  lERR  =  1 
RETURN 
END 
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